// 标准库或者第三方库

// 自定义头文件
#include "template_problem01.hxx"

// ACIS
#include "acis/acis.hxx"
#include "acis/api.err"
#include "acis/api.hxx"
#include "acis/body.hxx"
#include "acis/check.hxx"
#include "acis/cstrapi.hxx"
#include "acis/logical.h"
#include "acis/module.hxx"
#include "acis/primtive.hxx"

#include "acis/law_base.hxx"

#include<iostream>
#include<map>
#include<regex>
#include<vector>
#include<string>
#include<stack>
#include<complex>
#include<regex>
#include<cmath>
#include"../3rdparty/eigen/Eigen/Dense"
using namespace std;



/*
1. 首先判断是否为线性函数
2. 解析多项式字符串中的系数
3. 求解一元二次方程
4. 求解一元三次方程
5. 求解一元四-八阶方程 使用特征值求解 使用到Eigen库
6. 对三-八阶的结果进行带入验证，并使用牛顿迭代法降低误差
7. 
*/

// 存储多项式系数的数组
struct PolynomialResult {
	double* coefficients;  // 系数数组
	int size;             // 数组大小
	bool success;         // 转换是否成功
};

// 移除括号
static string removeBrackets(const std::string& input) {
	// 使用正则表达式替换所有左括号 '(' 和右括号 ')'
	string result = regex_replace(input, std::regex(R"([\(\)])"), "");
	return result;
}

// 解析系数 简单粗暴版本
static PolynomialResult parsePolynomialString(const string& polynomial_str) {
	PolynomialResult result;
	// TODO: 在这里实现解析逻辑
	result.coefficients = nullptr;
	result.size = 0;
	result.success = false;

	//int order;
	cout << "OriginData:" << polynomial_str << endl;
	string post_poly = removeBrackets(polynomial_str);
	// 最大阶数
	int max_order = post_poly[2] - '0';
	result.size = max_order + 1;
	result.coefficients = new double[max_order + 1](); // 全部初始化为0
	result.coefficients[0] = 1; // 第一个系数为1
	vector<size_t> x_pos;
	for (int j = 0; j < post_poly.size(); j++) {
		if (post_poly[j] == 'X'){
			// 存储每个x的位置
			x_pos.push_back(j);
		}
	}
	
	for (int i = 0; i < x_pos.size()-1; i++) {
		size_t coe_length = x_pos[i + 1] - x_pos[i] -3; // 这里有一个乘号
		//cout << "XXX" << coe_length << endl;

		result.coefficients[i + 1] = stod(post_poly.substr(x_pos[i]+3, coe_length));
	}
	// 常数系数
	result.coefficients[result.size - 1] = stod(post_poly.substr(x_pos[x_pos.size() - 1] + 1));

	//cout <<"First:" << result.coefficients[0] << endl;
	/*for (int k = 0; k < result.size; k++) {
		cout << "Coe:" << result.coefficients[k] << endl;
	}*/

	//cout << "Max Order:" << result.size << endl;


	return result;
}


//求解一元二次方程
static std::vector<double> solveQuadraticReal(const double* coef) {
	std::vector<double> realRoots;
	double a = coef[0], b = coef[1], c = coef[2];

	if (std::abs(a) < 1e-10) return realRoots;

	if (std::abs(a) < std::numeric_limits<double>::epsilon()) {
		// a≈0, 处理一次方程 bx + c = 0
		if (std::abs(b) > std::numeric_limits<double>::epsilon()) {
			realRoots.push_back(-c / b);
		} // else:  b也接近0，无解或无数解，返回空向量
		return realRoots;
	}

	double delta = b * b - 4 * a * c;

	if (std::abs(delta) < 1e-10) {  // 重根
		realRoots.push_back(-b / (2 * a));
	}
	else if (delta > 0) {  // 两个实根
		double sqrtDelta = std::sqrt(delta);
		realRoots.push_back((-b + sqrtDelta) / (2 * a));
		realRoots.push_back((-b - sqrtDelta) / (2 * a));
	} // else: delta < 0, 没有实根，返回空向量

	return realRoots;
}



class CubicEquation {
private:
	static constexpr double EPS = 1E-12;  // 计算精度
	static constexpr double SAME_ROOT_EPS = 1E-3;  // 重根判断精度
	static constexpr double SEARCH_RANGE = 100.0;  // 二分查找范围
	static constexpr double PI = 3.14159265358979323846;

	// 判断是否为0
	static bool isZero(double x) {
		return std::abs(x) < EPS;
	}

	// 计算函数值: ax³ + bx² + cx + d
	static double f(const double* coef, double x) {
		return ((coef[0] * x + coef[1]) * x + coef[2]) * x + coef[3];
	}

	// 计算导数值: 3ax² + 2bx + c
	static double df(const double* coef, double x) {
		return (3.0 * coef[0] * x + 2.0 * coef[1]) * x + coef[2];
	}

	// 二分查找
	static double binarySearch(const double* coef, double left, double right) {
		for (int i = 0; i < 100 && right - left > EPS; i++) {
			double mid = (left + right) / 2;
			double fmid = f(coef, mid);

			if (isZero(fmid)) return mid;

			if (f(coef, left) * fmid < 0) {
				right = mid;
			}
			else {
				left = mid;
			}
		}
		return (left + right) / 2;
	}

	// 牛顿迭代法精确化
	static double newtonRefine(const double* coef, double x) {
		for (int i = 0; i < 20; i++) {
			double fx = f(coef, x);
			if (isZero(fx)) break;

			double dfx = df(coef, x);
			if (isZero(dfx)) break;

			double dx = fx / dfx;
			x -= dx;

			if (std::abs(dx) < EPS) break;
		}
		return x;
	}

	// 求解二次方程
	static std::vector<double> solveQuadratic(double a, double b, double c) {
		std::vector<double> roots;

		if (isZero(a)) {
			if (!isZero(b)) {
				roots.push_back(-c / b);
			}
			return roots;
		}

		double discriminant = b * b - 4 * a * c;

		if (discriminant > 0) {
			double q = -0.5 * (b + std::copysign(std::sqrt(discriminant), b));
			roots.push_back(q / a);
			roots.push_back(c / q);
		}
		else if (isZero(discriminant)) {
			roots.push_back(-b / (2 * a));
		}

		return roots;
	}

	// 求解三次方程的解析解
	static std::vector<double> solveAnalytical(const double* coef) {
		std::vector<double> roots;

		double a = coef[0];
		double b = coef[1];
		double c = coef[2];
		double d = coef[3];

		// 标准化为 x³ + px + q = 0 的形式
		double p = (3.0 * a * c - b * b) / (3.0 * a * a);
		double q = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d) / (27.0 * a * a * a);

		// 计算判别式
		double D = q * q / 4.0 + p * p * p / 27.0;

		if (abs(D)<1e-10) {
			if (abs(p)< 1e-10) {
				// 三重根
				roots.push_back(-b / (3.0 * a));
			}
			else {
				// 一个单根，一个二重根
				double r = -3.0 * q / (2.0 * p);
				roots.push_back(-b / (3.0 * a) + r);
				roots.push_back(-b / (3.0 * a) - r / 2.0);
			}
		}
		else if (D < 0) {
			// 三个不相等的实根
			double theta = std::acos(-q / (2.0 * std::sqrt(-p * p * p / 27.0)));
			double r = 2.0 * std::sqrt(-p / 3.0);
			for (int k = 0; k < 3; k++) {
				roots.push_back(-b / (3.0 * a) + r * std::cos((theta + 2.0 * k * PI) / 3.0));
			}
		}
		else {
			// 一个实根，两个共轭复根
			double u = std::cbrt(-q / 2.0 + std::sqrt(D));
			double v = -p / (3.0 * u);
			roots.push_back(-b / (3.0 * a) + u + v);
		}

		return roots;
	}

public:
	// 主求解函数
	static std::vector<double> solve(const double* coef) {
		std::vector<double> roots;

		// 处理退化情况
		if (isZero(coef[0])) {
			return solveQuadratic(coef[1], coef[2], coef[3]);
		}

		// 解析解
		auto analytical_roots = solveAnalytical(coef);

		// 二分查找
		std::vector<double> binary_roots;
		double prev_f = f(coef, -SEARCH_RANGE);

		for (double x = -SEARCH_RANGE; x <= SEARCH_RANGE; x += 0.2) {
			double curr_f = f(coef, x);
			if (isZero(curr_f)) {
				binary_roots.push_back(x);
			}
			else if (prev_f * curr_f < 0) {
				binary_roots.push_back(binarySearch(coef, x - 0.2, x));
			}
			prev_f = curr_f;
		}

		// 合并所有找到的根
		roots.insert(roots.end(), analytical_roots.begin(), analytical_roots.end());
		roots.insert(roots.end(), binary_roots.begin(), binary_roots.end());

		// 对所有根进行精确化
		for (double& root : roots) {
			root = newtonRefine(coef, root);
		}

		// 排序并去重
		std::sort(roots.begin(), roots.end());
		auto last = std::unique(roots.begin(), roots.end(),
			[](double a, double b) { return std::abs(a - b) < SAME_ROOT_EPS; });
		roots.erase(last, roots.end());

		return roots;
	}
};




// 计算多项式值
double polynomial(const double* coef, int degree, double x) {
	double result = 0;
	for (int i = 0; i <= degree; i++) {
		result += coef[i] * pow(x, degree - i);
	}
	return result;
}

// 计算多项式导数值
double derivative(const double* coef, int degree, double x) {
	double result = 0;
	for (int i = 0; i < degree; i++) {
		result += (degree - i) * coef[i] * pow(x, degree - i - 1);
	}
	return result;
}

// 牛顿迭代法
double newton_method(const double* coef, int degree, double x0) {
	double x = x0;
	double prev_x;
	int iterations = 0;
	const int MAX_ITERATIONS = 10000;  // 最大迭代次数
	const double EPSILON = 1e-7;     // 误差阈值

	do {
		prev_x = x;
		double f = polynomial(coef, degree, x);
		double f_prime = derivative(coef, degree, x);

		// 防止除以0
		if (fabs(f_prime) < EPSILON) {
			//cout << "导数为0,迭代失败" << endl;
			return x;
		}

		x = x - f / f_prime;
		iterations++;
		if (iterations >= MAX_ITERATIONS) {
			return x;
		}
	} while (fabs(x - prev_x) > EPSILON);
	return x;
}

// 验证并优化根
void verifyAndOptimizeRoots(const double* coef, int degree, std::vector<double>& roots) {
	const double ERROR_THRESHOLD = 1e-6;

	for (size_t i = 0; i < roots.size(); i++) {
		double error = std::abs(polynomial(coef, degree, roots[i]));

		if (error > ERROR_THRESHOLD) {
			// 使用牛顿法优化此根
			//cout <<"优化前" << roots[i];
			roots[i] = newton_method(coef, degree, roots[i]);
			//cout << "优化后" << roots[i] << endl;
		}
	}
}


 //求解高阶次方程的根
class PolynomialSolver {
public:
	// 求解多项式方程，输入系数从高阶到低阶
	static std::vector<double> solve(const double* coef, int degree) {
		if (!coef || degree < 1) return std::vector<double>();

		// 构建伴随矩阵
		Eigen::MatrixXd companion = buildCompanionMatrix(coef, degree);

		// 计算特征值
		Eigen::EigenSolver<Eigen::MatrixXd> solver(companion);

		// 提取实特征值
		std::vector<double> roots;
		for (int i = 0; i < degree; i++) {
			std::complex<double> eigenvalue = solver.eigenvalues()(i);
			// 只保留实数根（虚部接近于0的特征值）
			/*if (std::abs(eigenvalue.imag()) < 1e-10) {
				roots.push_back(eigenvalue.real());
			}*/
			roots.push_back(eigenvalue.real());
			
		}
		/*for (int j = 0; j < roots.size() - 1; j++) {
			bool flag = true;
			for (int k = j + 1; k < roots.size(); k++) {
				if (std::abs(roots[j] - roots[k]) <= 1e-10) {
					flag = false;
				}
			}
			if (flag) {
				r_roots.push_back(roots[j]);
			}
		}*/

		//verifyAndOptimizeRoots(coef, degree, roots);
		if (!roots.empty()) {
			// 1. 先排序
			std::sort(roots.begin(), roots.end());

			// 2. 移除相邻的重复元素
			auto new_end = std::unique(roots.begin(), roots.end(),
				[](const double& a, const double& b) {
					return std::abs(a - b) < 1e-8;
				});

			// 3. 删除多余的元素
			roots.erase(new_end, roots.end());
		}

		return roots;
	}

private:
	// 构建伴随矩阵
	static Eigen::MatrixXd buildCompanionMatrix(const double* coef, int degree) {
		Eigen::MatrixXd companion = Eigen::MatrixXd::Zero(degree, degree);

		// 标准化系数（使最高次项系数为1）
		double leading_coef = coef[0];

		// 填充伴随矩阵的最后一行
		for (int i = 0; i < degree; i++) {
			companion(degree - 1, i) = -coef[degree - i] / leading_coef;
		}

		// 填充次对角线为1
		for (int i = 1; i < degree; i++) {
			companion(i - 1, i) = 1.0;
		}

		return companion;
	}
};



// ************************************************************
outcome answer_api_nroots_of_law(law* input_law, double start, double end, int* size, double** answer) {
	API_BEGIN
		if (api_checking_on) {
			check_law(input_law);
		}
    // 多项式字符串
    string temp_str = input_law->string();
	std::cout << input_law->string() << std::endl;
	//int dim = input_law->precedence();
	//cout << "DIM:" << dim << endl;

    int linear_flag = input_law->linear();
    if (linear_flag) {
		int a_pos=0;
		int b_pos=0;
		double a, b;
		int str_len = temp_str.length();
		 //首先判断第一个符号是
		for (int i = 0; i < str_len; i++) {
			if (temp_str[i] == 'X') {
				a_pos = i;
			}
			if (temp_str.substr(i,3) == "X^1") {
				b_pos = i+3;
			}
			
			break;
		}
		if (a_pos == 0) {
			a = 1;
		}
		else {
			a = stod(temp_str.substr(0, a_pos));
		}
		
		b = stod(temp_str.substr(b_pos));
		double result = - b / a;
		if (result >= start && result <= end) {
			//cout << "Result:" << result << endl;
			*size = 1;
			*answer = new double[1];
			(*answer)[0] = result;
		}


    }
	else {
		// 获取系数
		PolynomialResult result = parsePolynomialString(temp_str);

		// 4-8阶的方程求解
		if (result.size >= 5) {
			// 构造伴随矩阵
			//Eigen::MatrixXd companionMatrix = constructCompanionMatrix(result.coefficients);

			// 用 QR 分解提取特征值
			//vector<double> roots = solvePolynomial(result.coefficients, result.size-1);
			vector<double> roots = PolynomialSolver::solve(result.coefficients, result.size - 1);
			sort(roots.begin(), roots.end());
			*answer = new double[roots.size()];
			for (double root : roots) {
				if (root >= start && root <= end) {
					//cout << "root22:" << root << endl;
					(*answer)[*size] = root;
					*size = *size + 1;
				}

			}
		}

		// 求解一元三次方程
		if (result.size == 4) {
			*size = 0;
			std::vector<double> roots = CubicEquation::solve(result.coefficients);
			sort(roots.begin(), roots.end());
			*answer = new double[roots.size()];
			//std::cout << "实根数量: " << roots.size() << std::endl;
			for (double root : roots) {
				if (root >= start && root <= end) {
					(*answer)[*size] = root;
					*size = *size + 1;
				}
				
			}
		}

		// 求解一元二次方程
		if (result.size == 3) {
			*size = 0;
			//*answer = new double[result.size - 1];
			auto roots = solveQuadraticReal(result.coefficients);
			int count = 0;
			if (roots.size() > 0) {
				for (int i = 0; i < roots.size(); i++) {
					if (start < roots[i] && roots[i] < end) {
						count++;
						*size = *size + 1;
					}
				}
				// 从小到大排序
				sort(roots.begin(), roots.end());
				
			}
			if (count != 0) {
				*answer = new double[count];
				if (roots.size() > 0) {
					for (int i = 0; i < roots.size(); i++) {
						if (start < roots[i] && roots[i] < end) {
							if (count == 1) {
								(*answer)[0] = roots[i];
							}
							else {
								(*answer)[i] = roots[i];
							}
						}
					}

				}
			}
			
			
			
		}
	}


	API_END
		return result;
}
